Influence of counter-rotating von Karman flow 
on cylindrical Rayleigh-Benard convection 
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The axisymmetric flow in an aspect-ratio-one cylinder whose upper and lower bounding disks are 
maintained at different temperatures and rotate at equal and opposite velocities is investigated. In 
this combined Rayleigh-Benard/von Karman problem, the imposed temperature gradient is mea- 
sured by the Rayleigh number Ra and the angular velocity by the Reynolds number Re. Although 
fluid motion is present as soon as Re 7^ 0, a symmetry-breaking transition analogous to the onset 
of convection takes place at a finite Rayleigh number higher than that for Re — 0. For Re < 95, 
the transition is a pitchfork bifurcation to a pair of steady states, while for Re > 95, it is a Hopf 
bifurcation to a limit cycle. The steady states and limit cycle are connected via a pair of SNIPER 
bifurcations except very near the Takens-Bogdanov codimension-two point, where the scenario in- 
cludes global bifurcations. Detailed phase portraits and bifurcation diagrams are presented, as well 
as the evolution of the leading part of the spectrum, over the parameter ranges < Re < 120 and 
< Ra < 30 000. 

PACS numbers: 47.20.Ky, 47.20.Bp, 47.10.Fg, 47.32.Ef 



I. INTRODUCTION 

Thermal convection and shear are central to the study 
of hydrodynamic instabilities. The interest in these stems 
both from their large number of practical applications 
and also from their status as prototypes in theoretical 
investigation. In this study, we combine two well-known 
cylindrical configurations: Rayleigh-Benard convection 
in which the upper and lower bounding disks of the cylin- 
der are maintained at different temperatures, and von 
Karman flow in which the disks rotate at equal and oppo- 
site velocities. We choose the simplest small-aspect-ratio 
geometry: a cylinder with equal height and radius, with 
imposed axisymmetry. 

There is an extensive literature on Rayleigh-Benard 
convection in a small to medium aspect-ratio cylinder (in 
which the radius is between 0.5 and 5 times the height), 
e.g. [l]-!!)]. The references most relevant to us are axisym- 
metric simulations 7- 9] which produced patterns of ra- 
dially travelling concentric rolls (target patterns). These 
were shown to be produced in most cases by a saddle- 
node infinite period (SNIPER) bifurcation which will be 
the subject of much of the present investigation. 

The literature on von Karman flow is not as volumi- 
nous, but is growing. Restricting ourselves to the small 
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to medium aspect-ratio and the exactly or nearly exactly 
counter-rotating configuration, recent numerical articles 
include axisymmetric and non-axisymmetric simulations 
EMU- Interest in this problem has been enhanced by 
experiments using this configuration at high Reynolds 
numbers to generate turbulen ce fvA Il8j or a magnetic 
field via the dynamo effect |L9l - l21] . 

The combined Rayleigh-Benard/von Karman configu- 
ration seems, however, unexplored. This situation stands 
in contrast to that of other shear and rotating flows. 
The superposition of plane Poiseuille flow with Rayleigh- 
Benard convection (PRB), has been the subject of a 
great many investigations; indeed, a recent review [22| 
cites hundreds of articles. For longitudinal rolls in plane 
Poiseuille and Couette flows, the threshold remains un- 
changed [23j, |24j. For rolls transverse to the Poiseuille 
flow, the threshold increases and its nature changes from 
absolute to convective, e.g. [25|. In Taylor- Couette flow 
between cylinders at different temperatures, the nature 
of the instability also depends strongly on the angle be- 
tween the imposed shear and buoyancy force, e.g. (2a . l27j . 
Turning to rotating disks, heat transfer in the rotor- 
stator configuration (between a rotating and a station- 
ary disk) is extensively studied because of its applications 
to turbomachinery, e.g. [28]. Rotating Rayleigh-Benard 
convection has long interested geophysicists, as well as 
researchers in pattern formation e.g. [29j, because of the 
possibility of chaos at onset; SNIPER bifurcations are 
found in this system as well [3(| IHJ . Of the studies we 
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have found of the flow between counter-rotating disks 
maintained at different temperatures, one computes the 
similarity solution between disks of infinite radius prior 
to the onset of convection [32|, and another considers 
turbulent flows and neglects buoyancy (33j . 

The Rayleigh-Benard/von Karman problem in this ge- 
ometry turns out to display a number of fascinating phe- 
nomena. For this reason, this problem has been used to 
develop and test a reduced model for numerical simula- 
tion, in which fields are expanded in terms of the eigen- 
functions of the governing equations linearized about a 
steady state; this study is the subject of a companion 
paper [34j |. 

Our paper is organized as follows. Section ITTI presents 
the governing equations of our configuration, along with 
the non-dimensionalizations we have used. In section Hill 
we describe the various codes we have used to carry out 
time-dependent simulations, branch continuation and lin- 
ear stability analysis. Section IIVI shows how the transi- 
tion threshold is affected by the disk counter-rotation. In 
sections IV1 and fVTT we present the steady states and limit 
cycles resulting from the transition. The connection be- 
tween the steady states and limit cycle via a SNIPER 
bifurcation is explained in section IVlIl which presents a 
complete phase diagram and bifurcation analysis. Re- 
versing the traditional order, section IVIIII gives a linear 
stability analysis of the basic state, showing how the con- 
vective eigenvalues and eigenvectors are interlaced and 
joined by the von Karman flow. 



II. PROBLEM FORMULATION 



and lower bounding disks are thermally conducting, ro- 
tate with angular velocities +CI and —CI and are held at 
temperatures To and To + AT, respectively. The nondi- 
mensional parameters are defined as: 
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where v, a, 7 and g are the kinematic viscosity, the ther- 
mal diffusivity, the thermal expansion coefficient, and 
the gravitational acceleration, respectively. The Prandtl 
number Pr and aspect ratio T are set to one. Reynolds 
numbers in the range < Re < 120 will be studied, along 
with Rayleigh numbers in the range < Ra < 30000. 

Our study uses several independent codes, sum- 
marized in section Mil which use different non- 
dimensionalizations, and different problem formulations, 
i.e. primitive variables or streamf unction- vorticity. To 
facilitate later discussion, we briefly present these vari- 
ous formulations here. The non-dimensionalized temper- 
ature is taken to be the deviation from To divided by AT 
and lengths are non-dimensionalized by R. The imposed 
counter-rotation of the disks provides a natural timescale: 
if time is non-dimensionalized by I/O and velocities by 
RCl, the resulting governing equations are: 



v-u = 
<9 t U+(U-V)U = -VP 
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where A and V 2 denote the vector and scalar Laplacian, 
respectively. The boundary conditions are: 
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U r = U z = 0, Ug = +r, T = at z = l (3a) 

U r = U z = 0, U 6 = -r, T = 1 at z = (3b) 

U T = U Z = 0, U e = 0, d r T = at r = 1 (3c) 

U r = d r U z =0, U e = 0, d r T = at r = (3d) 

These equations cannot be used in the absence of ro- 
tation, when Re = 0. To allow for this case, the unit 
of time can instead be taken to be the viscous diffusion 
time R? /v, with the resulting velocity scale v/R. In ad- 
dition, with the assumption of axisymmetry, the merid- 
ional velocity components (V r , V z ) are best described by 
the Stokes streamf unction VP: 



FIG. 1: Cylinder with height and radius equal to one. The 
bottom disk is heated and rotated in the clockwise direc- 
tion, while the top disk is cooled and rotated in the counter- 
clockwise direction. 

We consider a viscous Newtonian fluid governed by the 
Boussinesq approximation and contained in a cylinder of 
radius R and height H. The velocity and temperature 
fields are assumed axisymmctric. The bounding cylin- 
der is stationary and thermally insulating. The upper 
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the azimuthal vorticity is: 



x V* 



T> 2 * = - 



- -cO * = - (d-d r 
r I r 



(4a) 



(4b) 



+ dl)^ (4c) 



3 



and the Navier-Stokes and Boussinesq equations reduce 
to: 

d s D 2 V + (V r d- + V z d z ) D 2 ^ = ( V 2 - 1) L> 2 * 



—d r T + d z M- 
Pr \ r 

1 



(5a) 



d s T+(V r d r + V z d z )T = — V 2 T (5b) 

Pr 

d s V e + (V r d+ + V z d z ) V = (V - V e (5c) 

The rotation then appears in the boundary conditions, 
which become: 

* = d z y = o <=> v r = v z = o, 

Vg = +rRe, T = at z = 1 (6a) 

V = d x V = 0<=>V r = Vz = 0, 

Vg = -rRe, T = 1 at z = (6b) 

V = d r V = 0<=>Vr = V z = 0, 

Vg = 0, <9 r T = at r = 1 (6c) 
* = L> 2 * = K r = d r V z = 0, 

= 0, d r T = at r = (6d) 

The velocity and time in (J2j)-(j3j) are related to those 
in ©-© by V = i?eU and s = t/Re; the temper- 
ature T and all lengths remain unchanged. Equations 
([5c]) and (|5bl) can be considered to be advection-diffusion 
equations for the azimuthal velocity Vg and the tem- 
perature T. Equation (|5a|) shows that gradients in Vg 
and T in turn generate vorticity. The vertical gradient 
d z (Vg 2 /r), when combined with the boundary conditions 
(|6ap - (|6b[) . corresponds to Ekman pumping, by which the 
z-dependent Vg generates meridional vorticity and veloc- 
ity for any non-zero Re. This meridional velocity in turn 
affects the azimuthal velocity and the temperature via 
equations (|5c|) and (I5bp ." However, the thermal gradi- 
ent term d r T in the vorticity equation (|5ap is radial, and 
so the axial temperature gradient imposed by boundary 
conditions in (|6ap ~ (l6bp docs not immediately impart mo- 
tion to the fluid for any non-zero Ra. That is, Rayleigh- 
Benard convection occurs via an instability at a finite 
threshold Ra c , in contrast to natural convection, in which 
an externally imposed horizontal gradient generates mo- 
tion for any Ra > via the term Ra d r T. 

Finally, we state the symmetries of the configuration 
of our problem, which can be seen from figure [T] The 
Rayleigh-Benard convection problem has the symmetry 
of combined reflection in z and in T about their mean 
values, sometimes called the Boussinesq symmetry. The 
von Karman flow has the symmetry of combined reflec- 
tion in z and in 9, called R^ [15j. (In this axisymmetric 
problem, reflection in 9 means only reversing the sign of 
Ug; since all solutions are invariant under rotations in 
9, these are not considered.) The Rayleigh-Benard/von 
Karman problem thus has the symmetry which combines 
reflection in z, 9 and T . A state is reflection-symmetric 
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if it is invariant under the reflection operator defined in 
the two formulations by: 



(r, 1 — z) or 



n:\Ug] (r,z)= -Ug | (r, 1 - z) (7) 
T / V 1-T 



III. NUMERICAL METHODS 

We have solved the governing equations using a num- 
ber of different codes. The first is the time-integration 
code written by Patankar [35[ for equations (HJ)-©. The 
temporal discretization uses the backwards differentia- 
tion formula for the time derivative, the explicit Adams- 
Bashforth formula for the advective and buoyancy terms, 
and implicit evaluation of the diffusive terms, which leads 
to the following form for timestepping equation @: 
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The spatial discretization uses centered finite differences. 
A staggered grid is used, with pressure and temperature 
represented on the nodes and the velocity components 
between the nodes. A timestep At = 1CP 3 and a reso- 
lution in the radial and axial direction of 82 x 82 were 
used. 

Time integration has also been carried out with a sec- 
ond code, which uses the streamfunction-vorticity for- 
mulation ([S])-®, a second-order finite difference scheme 
in space, and an alternating-direction-implicit scheme in 
time [36(. A 51 x 51 grid with a typical timestep of 
As = 10 -4 were found to be adequate for obtaining con- 
verged results. 

We also carried out bifurcation and linear stability 
analysis. We followed branches via Newton's method 
and calculated eigenvalues via ARPACK. Several differ- 
ent spatial representations have been used. One code uses 
the same spatial discretization - i.e. the streamfunction- 
vorticity formulation and the same spatial resolution - 
as the time-integration code mentioned in the previous 
paragraph. The Jacobian of the matrix built in the last 
step of Newton iteration is directly used in ARPACK in 
shift-invert mode. The branches are followed using an 
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arclength continuation as described in [37| and the prac- 
tical implementation of bordered matrix inversion can be 
found in 

The second code (39j with which we carried out bi- 
furcation and linear stability analysis uses the primitive 
variable formulation and a pseudo-spectral spatial dis- 
cretization. The number of Chebyshev polynomials (or 
Gauss-Lobatto collocation points) used is 15 in each of 
the r and z directions. 

The quantitative results which we present have been 
obtained by two or more of the codes described above, 
and have also been verified by varying the resolution. 



IV. THRESHOLDS 

30000 



20000 - J- 




o - 



I I I I I I I I I I I I I I I I I I I I I I I 

20 40 60 80 100 120 

Reynolds 

FIG. 2: (Color online) Convection thresholds as a function of 
Reynolds number. Solid blue curve: transition to steady con- 
vection. Dashed red curve: transition to oscillatory convec- 
tion. Dotted black curve: fit to even fourth-order polynomial 
@. Points: thresholds obtained from time-dependent code. 
Letters indicate parameter values for flows presented in fig- 
ures [3] and [4] (B, basic state), [5] and [6] (S, steady convection), 
and [7] and [5] (O, oscillatory convection). 

In the absence of rotation, i.e. Re = 0, the govern- 
ing equations have as a solution the conductive state, in 
which the fluid is stationary and the temperature varies 
linearly in the vertical direction: 

U = 0, T = 1 — z (8) 

Although the conductive solution exists for all Ra, it loses 
stability at some critical value Ra c to a convective solu- 
tion, in which fluid motion ensues, breaking the Boussi- 
nesq reflection symmetry (J7J . In the von Karman system, 



as mentioned in section [TT1 Ekman pumping leads to re- 
circulating regions for any non-zero value of Re. How- 
ever, for a fixed value of Re, there is still a transition at 
a well-defined Ra c , at which one solution loses stability 
to another. By analogy with the non-rotating case, we 
will call these the basic and convective solutions, despite 
the presence of fluid motion (and thus of convective heat 
transfer) for Ra < Ra c , and we will continue to refer to 
the transition between them as the onset of convection. 

We begin by showing the critical Rayleigh number 
for onset of axisymmetric convection as a function of 
Reynolds number. For Re < 95, the bifurcation is to 
a steady state (i.e. the critical eigenvalue is real), and 
for Re > 95, the bifurcation is to an oscillatory state 
(i.e. the critical eigenvalue is complex). The thresholds 
in figure [2] are obtained by calculating the eigenvalues 
of the linearized problem for fixed values of (Re,Ra) 
and interpolating to find the value Ra c (Re) at which the 
eigenvalue or its real part crosses zero. Thresholds were 
also obtained independently by fitting decay rates from 
the time-stepping code to exponentials and extrapolating 
these to zero. 

Away from where the bifurcation changes its nature 
from steady to oscillatory, Ra c should be a smooth func- 
tion of Re. In fact, Ra c must be an even function of Re, 
since reversing Re means reversing the direction of ro- 
tation of the upper and lower disks, an operation which 
leaves all dynamical properties unchanged, as was argued 
[2o| in the context of Rayleigh-Benard convection with 
throughflow. We fit Ra c over the range of the steady 
bifurcation, < Re < 94, to an even polynomial, obtain- 
ing: 

Ra c « 2260 + (0.6243 x Re) 2 + (0.0840 x Re) 4 (9) 

We see that the von Karman flow stabilizes the system 
against thermal convection. We contrast this with re- 
sults from other mixed thermal/shear systems. In plane 
Poiseuille and Couette flow, the threshold for onset of 
longitudinal rolls is unaffected by the shear [23|, [24[ . (In 
our system, the rolls can be seen as longitudinal, since the 
axes of concentric rolls are azimuthal, as is the velocity of 
the bounding disks.) In Taylor-Couette flow (measured 
by the Taylor number Ta) with a radial temperature gra- 
dient (measured by the Grashof number Gr), the role of 
rotation and heating are reversed from our system. Any 
radial temperature gradient causes large-scale motion (a 
role played by Re in our system) and a bifurcation to 
Taylor vortices is seen at a finite rotation rate (a role 
played by Ra in our system). In experiments on this 
system [26j, Ta c decreases with Gr, in contrast to the 
increase in Ra c with Re in our case. 

The threshold value Ra c (Re = 0) = 2260 agrees well 
with those previously reported for the onset of convec- 
tion for insulating sidewalls and aspect ratio F = 1, 
notably 2260 [J, Q, 2300 (ARa c /Ra c =1.8%) 4], 2241 
(ARa c /Ra c = 0.8%) @, and 2250 (ARa c /Ra c = 0.4%) 
Q. These authors all found the most unstable mode to 
be axisymmetric for cylindrical convection at this aspect 



FIG. 3: (Color online) Basic flow at Re = 40, Ra = 2000. 
From left to right: Ug, 9 and T. Ranges for U r and U z are 
[-0.05,0.05] and [-0.04,0.04], respectively 



FIG. 5: (Color online) Convective state at Re = 40, Ra = 
4000. From left to right: Ug, * and T. Ranges for U r and U z 
are [-0.11,0.19] and [-0.29,0.11], respectively. 




FIG. 4: (Color online) Basic flow at Re = 90, Ra = 8000. 
From left to right: Ue, 9 and T. Ranges for U r and U z are 
[-0.11,0.11] and [-0.07,0.07], respectively. 
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FIG. 6: (Color online) Convective state at Re = 90, Ra = 
12500. From left to right: Ug, * and T. Ranges for U r and 
U z are [-0.15,0.18] and [-0.24,0.10], respectively. 



ratio. 

Varying Re and investigating the thresholds of non- 
axisymmctric cigenmodes shows that the critical az- 
imuthal wavenumber is m = 2 for Re > 54.8. At Re = 0, 
the primary axisymmetric convective branch becomes un- 
stable to an to — 2 perturbation near Ra — 3000 [3, Q . 
For the pure von Karman flow, i.e. with Ra set to zero 
and Re gradually increased from zero, the first instabil- 



ity occurs near Re = 300 to an to = 2 mode [15j. We 
nonetheless restrict our consideration here to the axisym- 
metric bifurcations and eigenmodes. There are a num- 
ber of reasons for doing so. First, an understanding of 
the axisymmetric problem is important for completeness: 
branches which are unstable may nonetheless play a role 
in the dynamics. Secondly, the axisymmetric problem 
displays a number of fascinating phenomena, both lin- 
ear and nonlinear, as we will see. Because of this, it has 
served as a test case for a numerical method using a re- 
duced model [HI ; a presentation of the results from fully 
resolved computations is a necessary complement to this 
study. Finally, additional physical effects, e.g. magnetic 
fields, can stabilize the axisymmetric configuration (40| . 



V. STEADY CONVECTIVE STATES 

We now discuss the nature of the steady flows below 
and above the convection threshold. The basic state, 
which is the analogue of the conductive state in the pres- 
ence of von Karman flow, is shown in figures [3] and U for 
Re = 40, Ra = 2000 and for Re = 90, Ra = 8000 respec- 
tively. The azimuthal velocity increases gradually from 
negative in the lower half to positive in the upper half 
of the cylinder, following the counter-rotating disks. It is 



singular at the corners r = 1, z = 0, 1 since the boundary 
conditions are discontinuous where the stationary cylin- 
der meets the rotating disks. The large recirculating cells 
caused by Ekman pumping are prominent features in the 
(r, z) plane, carrying fluid radially outward along both ro- 
tating disks at z = 0, 1, then along the bounding cylinder, 
inwards along the midplane z = 1/2, and back towards 
the two disks along the axis. The azimuthal velocity 
and recirculating cells combine to yield fluid trajectories 
which are toroidal. 

We may compare the basic flows in figures [3] and 2] 
For Re = 90, the temperature field is noticeably dif- 
ferent from the linear conductive profile ([5]). The re- 
circulating cells are stronger, with a maximal velocity 
of |C/, max | = 0.11, as compared to |C/, max | = 0.04 for 
Re = 40. These stronger recirculating cells exert more in- 
fluence on the temperature field. At r — 1, the isotherms 
are deviated towards the midplane as warm fluid con- 
verges upwards from the lower disk and cold fluid down- 
wards from the upper disk. At r = 0, the isotherms 
are deviated again, here outwards towards the upper and 
lower disks. The azimuthal velocity gradients for Re = 90 
are more concentrated near the disks. The basic states 
shown in figures [3] and E] are reflection-symmetric about 
the midplane: each is invariant under defined in 0. 

The transition to convection breaks the reflection sym- 
metry. In the absence of rotation, the transition occurs at 
Ra = 2260, as stated in section ITVl and leads to a single 
roll which occupies the entire cavity. For Re = 40, tran- 
sition occurs at Ra c = 3 125 and the resulting convective 
state, shown in figure [5] for Ra — 4000, has one small roll 
(due to the weak von Karman flow) in the upper right 
corner in addition to the main large convection roll. For 
Re — 90, the onset of convection occurs at Ra = 8 622 
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t = t = 4.5 ~ t/S t = 10.5 t = 19.4 t = 29.5 ~ 3r/4 t = 33sT 

FIG. 7: Near-sinusoidal limit cycle at Ra — 20 000, Re = 110. \E r is shown at times near 0, t/S, t/4, r/2, 3r/4, r indicated on 
timeseries in figure[9] Times are given in units of I/O.. The overall limit cycle has reflection symmetry: the states in the second 
half of the cycle are related by reflection symmetry to those in the first half of the cycle. The vortex in the upper right corner 
is smaller during the first half of the cycle and larger during the second half. The streamfunction contours are not equally 
spaced, but instead chosen to illustrate topological features of the flow. 




j. j, j, j. j, T 



t = 23.a = *d t = 68.7 t = 71.0 t = 72.0 t = 73.0 t = 86.5 s id + r/2 

FIG. 8: Near-heteroclinic limit cycle at Ra — 20 000, Re = 63. ^ is shown at times indicated on timeseries in figure^ t — 22.6; 
four instants in 68 < t < 73, during which the flow changes a great deal; and t = 86.5, approximately a half-period r/2 « 57.5 
after t = 22.6. The vortex in the lower right corner grows and the vortex in the upper left corner folds and then divides into 
three small vortices. Two of these disappear, leaving a small vortex in the upper right corner. As in figure [7] the states in 
the second half of the cycle are related by reflection symmetry to those in the first half of the cycle and the streamfunction 
contours are chosen to illustrate topological features of the flow. 



and the resulting convective state is shown in figure [5] at 
Ra = 12 500. The influence of the counter-rotating disks 
is stronger than it is at Re = 40 and so the recirculat- 
ing rolls, resulting from the combined effects of Ekman 
pumping (which favors two equal rolls) and convection 
(which favors a single roll), are more equal in size. The 
consequences on Ug and T of these changes in 'F can also 
clearly be seen in figures \5\ and [6] (These in turn affect 
'F via the buoyancy and Ekman pumping effects, though 
to a lesser extent.) 

Although the smaller cell in figures [5] and [5] is in the 
upper corner of the cylinder, the pitchfork bifurcation can 
produce a small cell at either the upper or lower corner. 
The breaking of reflection symmetry is also manifested 
by the azimuthal velocity: the Ug = contour intersects 
the axis near z = 0.2 in both figures [5] and [6j in contrast 
to the basic flow, for which this contour is a straight line 
across the midplane at z = 0.5. 



VI. OSCILLATORY CONVECTIVE STATES 

For Re > 95, the transition from the basic state is a 
Hopf bifurcation, creating oscillatory flows. We fix Ra = 
20 000 and present two oscillatory states, or limit cycles, 
one at Re = 110 in figure [7] and one at Re = 63 in 



figure [51 (Although for Re = 63, the first bifurcation as 
Ra is increased is to steady convection, there is another 
bifurcation at higher Ra to oscillatory convection and the 
resulting limit cycles are stable and connected to those 
created by the Hopf bifurcation for Re > 95, as we will 
see in section IVII1 ) We call the limit cycles at Re — 
110 and Re = 63 near-sinusoidal and near-heteroclinic, 
respectively, which we will also explain in section IVIII 

Figures [7] and [5] show meridional vortices jostling with 
one another. At Re=110 (figure[7]), the vortices are some- 
what similar in size, with the vortex in the upper or lower 
outer corner alternatively becoming smaller and larger. 
At Re=63 (figure [5]) , the principle vortex remains much 
larger throughout the limit cycle, with much smaller vor- 
tices on the periphery. At t = 71.0, there appear to be 
two vortices which both extend over the entire height, 
suggesting an interpretation of this limit cycle as a com- 
petition between one and two concentric radial vortices. 
This interpretation will be discussed further in section 

Enn 

Section [V] showed that the transition to steady con- 
vection breaks the reflection symmetry, creating two 
asymmetric states related by the reflection operator k. 
The transition to oscillatory convection has an analogous 
property: each of the instantaneous states shown in fig- 
ures [7] and [5] is asymmetric, and the second half of each 
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20 40 60 80 100 120 140 0.2 0.4 0.6 0.8 

time T(r=0.16,z=0.5) 



FIG. 9: Timeseries (left) and phase portraits (right) corresponding to limit cycles at Ra = 20 000. Dots in timeseries refer 
to visualizations in figures [7] and [8] Dots in phase portraits are plotted at equally spaced times so that their density reflects 
the rate at which the limit cycle is traversed. Top: at Re = 110, the oscillations are regular and near-sinusoidal. Bottom: at 
Re = 63, the oscillations have a much longer period and consist mainly of two long plateaus with abrupt gradients between 
them. 
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FIG. 10: Variation of square frequency with Re. Left: 1/r 2 for Ra = 20 000 (hollow dots) and for Ra = 18 000 (solid dots). 
Right: (Re/r) 2 for Re = 20 000 (hollow dots) and for Ra = 18 000 (solid dots). 



limit cycle is related to the first half by reflection, i.e.: 

(U r ,U e ,U z ,T)(t + T/2) = K{U r ,U e ,U z ,T)(t) (10) 

where r is the oscillation period. 

To better understand the limit cycles in figures [7] 
and[5J we plot various scalar quantities - l/ z (0.16, 0.5), 
U e (0.08, 0.5) and T(0.16,0.5) - as a function of time and 
of each other to create the timeseries and phase portraits 
of figure Figure[S] clearly shows the very different char- 
acter of the limit cycles at Re = 110 and Re = 63. At 
Re = 110, the timeseries is fairly close to harmonic, and 
the phase portrait shows that the limit cycle is traversed 
at a fairly constant speed. In contrast, at Re — 63, the 
timeseries shows two phases during which change is very 
slow, punctuated by abrupt transitions. The phase por- 
trait corroborates this: two regions of the limit cycle show 



a dense accumulation of points, indicating a slow traver- 
sal of these regions. 



We have seen that the period at Re = 63 is about three 
times that at Re — 110; it in fact approaches infinity as 
Re is decreased. (The approach by r to infinity implies 
that it is extremely sensitive to any change in the physical 
or numerical parameters. Thus, values of r are subject 
to a great deal of uncertainty.) This is demonstrated 
in figure [TU1 which plots the dependence of the square 
frequency (inverse square period) on Re for Ra = 20 000 
and for Ra = 18 000, for two different scalings: with 
the disk rotation period l/fl (left) and with the viscous 
diffusion time R 2 jv (right). 
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VII. BIFURCATION DIAGRAM 
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FIG. 11: (Color online) Curves of bifurcation points in the 
(Re,Ra) plane showing first pitchfork (blue, solid, PFi), sec- 
ond pitchfork (green, short-dashed, PF2), Hopf (red, long- 
dashed, H), and saddle- node (black, dotted, SN) bifurca- 
tions. Dots indicate codimension-two points, where curve 
PF 2 meets curves SN (90,14300) and H (95,11750). An 
enlargement of the region inside the square is shown in fig- 
ure 1121 Phase portraits shown as insets, a, b) Below PFi 
and H, the only solution is the stable basic flow, which is a 
node (a) or a spiral focus (b). c) Between PFi and PF2, the 
basic flow is unstable and there exist two stable asymmetric 
convective states, d, e) Between PF2 and SN, there exist 
five states: the unstable basic flow and two stable and two 
unstable convective states, f, g) Between SN and H, the un- 
stable basic flow is surrounded by a stable limit cycle, which 
is near-heteroclinic near SN (f) and near-sinusoidal near H 
(g)- 

We have mentioned several features of the oscillatory 
states in section IVT1 

• The small vortex does not remain in the upper or lower 
corner, as it does for the steady convective states, but al- 
ternates between the two locations, as shown in figures [7J 
and [8] 

• The character of the oscillations varies substantially, 
from near-sinusoidal to relaxational (plateaus punctu- 
ated by fast changes) as shown in figure 

• The period varies substantially with changes in Re or 
Ra, as shown in figure [TUl 

All of these features are typical of the bifurcation sce- 
nario we will now describe, and which is illustrated in 
figure [TTJ 

For Re < 95, the threshold for convection, already 
shown in figure [H corresponds to a pitchfork bifurcation 



PF\. For Ra below PF±, i.e. inset (a) in figure [TTJ the 
only solution is the basic flow, an example of which is 
given in figure[3l Proceeding clockwise, at PFi, the basic 
state loses stability and gives rise to two symmetrically- 
related convective states, as illustrated in inset (c) and 
figure[5l A second pitchfork bifurcation, PF2, occurs at a 
higher value of Ra. For Ra between curves PF% and SN, 
in regions (d) and (e), there exist two additional steady 
states, both unstable and also related by the symmetry 
operation ([7]). The trajectories leaving the two unstable 
states terminate on the two stable states, forming an in- 
variant set that resembles a circle or ellipse. Leaving PF2 
and approaching SN, as shown in (e), the two stable and 
unstable states approach one other along this invariant 
set, finally merging and annihilating one another in two 
symmetrically related saddle-node bifurcations at SN. 

Exactly at SN, trajectories spend an infinite period of 
time at the two locations at which the stable and unsta- 
ble states merged, forming a limit cycle whose period is 
infinite: a heteroclinic cycle. For Ra above SN, but suf- 
ficiently nearby, as in inset (f), much of the time during 
the limit cycle is spent in the neighborhood of the for- 
mer fixed points, as shown for Re = 63, Ra = 20 000 in 
figure [5] and in the bottom half of figure [51 Leaving SN, 
by increasing Ra or Re, as in inset (g), the time spent 
in the vicinity of these points shortens, and the period 
decreases, as shown for Re = 110, Ra = 20 000 in figure 
[7] and in the upper half of figure [5J Finally, for Re > 95, 
curve H corresponds to a Hopf bifurcation, below which 
the basic state is a stable spiral focus, as in inset (b). 

Although the saddle-node bifurcation signals the merg- 
ing of two pairs of steady states, it results in a limit cycle 
because of the connections between these steady states. 
(These connections in turn result from the formation of 
the steady states via two successive pitchfork bifurca- 
tions.) Under these circumstances, the saddle-node bi- 
furcation is called a SNIPER, for Saddle-Node In a PERi- 
odic orbit, or Saddle-Node Infinite PERiod; other names 
are a saddle-node homoclinic, a SNIC (Saddle-Node on 
Invariant Circle), or an Andronov bifurcation [u], l42l |. 
It can be shown that, near a SNIPER bifurcation, the 
period t of the limit cycle varies like 

1 . . 

where fi — fi c can be either Ra — -RagN for fixed Re or 
Re — i?esN for fixed Ra or any combination of the two. 
Although the square frequency is proportional to the dis- 
tance from threshold sufficiently near the threshold for 
any scaling, figure [TUl shows that, for our case, the linear 
dependence (jlip holds over a wider range of Re when t 
is scaled by the viscous diffusion time (right) than when 
it is scaled by the disk rotation period (left). 

Figure [TTJ shows two codimension-two points. At 
(95,11750), curves H, PF\ and PF 2 meet in a Takens- 
Bogdanov (TB) point. At (90,14300), curve SN ter- 
minates on PF2 at a hysteresis point, where the pitch- 
fork bifurcation changes from supercritical to subcritical. 
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FIG. 12: (Color online) Left: enlargement of square in figure 1111 The additional curve (violet, dash-dotted) indicates the 
secondary Hopf bifurcation SH which exists between codimension-two points (95, 11 750) (where it meets H) and (74, 17 100) 
(where it meets SN). Labels correspond to those of figure ITT1 Behavior in the square is described by the TB normal form l|12fl. 
Right: behavior of the TB normal form (|12p . Shown are the pitchfork bifurcation P, Hopf bifurcation H, secondary Hopf 
bifurcation SH, as well as the saddle-node of periodic orbits SNP and the gluing bifurcation G which are contained in the 
unfolding of the TB point. Phase portraits b), c), g) as in figure QTJ h) Between P and SH a stable limit cycle surrounds the 
unstable basic and convective flows, i) Between SH and G, small unstable limit cycles surround each stable convective state, 
j) Between G and SNP, two large limit cycles, one stable and one unstable, surround the three steady states. Hatched region 
indicates bistability between limit cycle and steady states. 



(Because the TB point was calculating using eigenvalue 
computations rather than by specialized algorithms such 
as those in (38[, we have been able to determine its lo- 
cation only to within about 0.5% in Re and 5% in Ra.) 
Between these two points, figure QTJ shows curve PF 2 
separating region (c), in which the unstable basic state 
coexists with two stable steady states, and region (f,g) in 
which it coexists with a stable limit cycle. The pitchfork 
bifurcation PF 2 cannot bring about such a transition. 
(The other paths between (c) and (f,g), either clockwise 
through PF 2 , region d/e and SN; or counter-clockwise 
through PF\, region a/b and H, can bring about this 
transition.) This means that figure QTJ is incomplete. 
The solution is to be found in the normal form of the 
Takens-Bogdanov codimension-two point in the presence 
of reflection symmetry [42|, |43| : 



dx 

~db 
dy 

dt 



x 2 y 



(12a) 
(12b) 



where 6 = ±1. We take S — +1, since it is this choice 
which reproduces the phenomena we see in the fluid- 
dynamical simulation. 

The behavior of the normal form (|12p is illustrated in 



Steady state 


Jacobian 


Eigenvalues 


(0,0) 

(±V=mT,o) 


(° ') 

\ -Mi M2 / 
( 1 \ 
y 2/Ui M2M1 / 


A± = §M2 ± \Jjf4 - Mi 
A± = |(M2 + Mi) 





TABLE I: Properties of the Takens-Bogdanov normal form. 



figure Q2] (right) and in table|H The TB codimension-two 
point is located at (mi,M2) = (0,0) and the basic state 
at (x, y) = (0,0). The basic state undergoes a pitchfork 
bifurcation (A + = 0), analogous to our PF\ and PF 2 , 
at Mi = and a Hopf bifurcation (TZe(X±) — 0), analo- 
gous to our H at fi 2 — for mi > 0. Additionally, table 
U shows that the solutions (±\/— Mi> 0) undergo a Hopf 
bifurcation at fi 2 — —(Jti for Mi < 0. Indeed, we have 
calculated a secondary Hopf bifurcation SH for the hy- 
drodynamic system, and indicate it in figure 1121 This 
bifurcation is subcritical and the limit cycles it gener- 
ates are unstable. Arnold [iHEil has shown that system 
(|12p contains two additional bifurcation curves emanat- 
ing from the codimension-two point (pi, fi 2 ) = (0, 0), one 
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a gluing bifurcation, G, and the other a saddle-node of 
periodic orbits, SNP. Although the curves SNP and G 
correspond to global bifurcations and cannot be calcu- 
lated by examination of the eigenvalues in table [IJ they 
have been determined [42| to be [12 = —0.752 /ij for SNP 
and 1x2 = — 0-8 \i\ for G. (We have computed neither for 
our hydrodynamic system.) Together, bifurcations SNP, 
G, SH and P accomplish the transformation in the small 
square of figure [TTJ between the three steady states of re- 
gion (c) to the stable limit cycle of region (g) . 

We describe the bifurcations and phase portraits in fig- 
ure H21 proceeding around the TB point in the counter- 
clockwise direction. In region (b) , the only solution is the 
stable basic state. The Hopf bifurcation H destabilizes 
the basic state and generates a stable limit cycle (g) . The 
pitchfork bifurcation P is subcritical, and generates two 
unstable asymmetric states while reducing the number 
of unstable directions of the basic state from two to one 
(h). The secondary Hopf bifurcation SH is also subcriti- 
cal, stabilizing the two asymmetric states while creating 
unstable limit cycles encircling each one (i). The gluing 
bifurcation G joins these two limit cycles to the symmet- 
ric state in a figure-eight which, afterwards, form a single 
unstable limit cycle encircling all three steady states (j). 
Curves SH and SNP delineate a narrow wedge in which 
there is bistability between the limit cycle and the steady 
convective states. The stable and unstable limit cycles 
annihilate in a saddle-node of periodic orbits SNP, leav- 
ing only the unstable basic flow and the stable asymmet- 
ric steady states (c). In crossing the pitchfork bifurca- 
tion P, we return to the single stable symmetric state of 
b). Figure [T2l shows a third codimension-two point where 
SH and SN meet at (74,17 100). Other codimension-two 
points which we have not located must mark the ends of 
curves G and SNP as well. 

Although the analysis surrounding the TB point is im- 
portant for completeness, we emphasize that, over most 
of the (Re, Ra) plane, it is the SNIPER bifurcation which 
marks the boundary between steady and oscillatory con- 
vection; other bifurcations come into play only between 
(95,11750) and (74,17100). 

Figure [T3] shows representative bifurcation diagrams. 
Branches are labelled by their stability index, i.e. the 
number of eigenvalues with positive real parts. The sta- 
bility index changes by one at a pitchfork bifurcation (for 
the basic state) or a saddle-node bifurcation and by two 
for a Hopf bifurcation. We plot the temperature at a 
fixed point T = T(r = 1/2, z = 1/2) as a function of 
Ra for fixed Re along the bottom of the figure, and as a 
function of Re for fixed Ra along the right. To represent 
limit cycles, we plot [~ J r T(t) 2 dt] 1 / 2 , the L 2 norm of 
T over an oscillation period r. The three codimension- 
two points at (74, 17100), (90, 14300) and (95,11750) 
and the limits of the pitchfork and saddle-node curves 
(Ra = 2260 and 28 445) delimit the ranges over which 
different types of bifurcation diagrams occur. 
• Re = 60 represents the large typical range < Re < 74: 
two successive supercritical pitchfork bifurcations gener- 



ate two pairs of steady states which are annihilated by 
a saddle-node (SNIPER) bifurcation, leading to a limit 
cycle. 

• Re — 85 represents [74, 90]: secondary Hopf bifurca- 
tions have been added to the diagrams. 

• Re — 93 represents [90, 95]: the second pitchfork bi- 
furcation has become subcritical. 

• For Re > 95 (not shown), there are no non-trivial 
steady states and a Hopf bifurcation leads directly to 
the limit cycle. 

• Ra = 5000 represents the range 2260 < Ra < 6640: 
a single pitchfork bifurcation leads to a pair of steady 
states as Re decreases. 

• Ra = 7000 represents [6640, 11750]: an additional 
pitchfork bifurcation and pair of steady branches can be 
seen. 

• Ra = 13 000 represents [11 750, 14 300]: a primary and 
two secondary Hopf bifurcations have been added to the 
scenario. 

• Ra = 15 000 represents [14 300, 17 100]: the two pairs of 
non-trivial branches are no longer connected to the triv- 
ial branch via pitchfork bifurcations, but instead arise via 
saddle-node bifurcations. 

• Ra= 18 000 represents [17100, 28 445]: the secondary 
Hopf bifurcations have disappeared. 

• For Ra > 28 445 (not shown), there are no non-trivial 
steady states, but only a limit cycle generated by the 
Hopf bifurcation. 

(The other bifurcations depicted in the enlargement of 
figure [121 which are present for the ranges 74 < Re < 95 
and 11 750 < Ra < 17 100, are omitted from figure [TBI for 
clarity. ) 

Scenarios similar to those reported in this section have 
been computed in other hydrodynamic configurations. 
We review these studies and list the similarities and dif- 
ferences between these observations and the current re- 
sults. In a cylinder with aspect ratio T = 5, Pr = 10 
and thermally conducting sidewalls @, H| , steady states 
are produced by two pitchfork bifurcations. When Ra 
is increased, a SNIPER bifurcation leads to a pattern of 
five concentric rolls traveling radially inwards. By reduc- 
ing the sidewall conductivity, other saddle-node bifur- 
cations on the same invariant circle appeared, creating 
stable steady four-roll states. No Hopf bifurcations or 
Takens-Bogdanov points were found. Such a bifurcation 
may also have been responsible for long-period oscilla- 
tions observed in experimental timeseries, but for which 
the flow could not be visualized [2]. A SNIPER bifur- 
cation is also found for r = 4, Pr = 7 in a study of 
rotating convection [30j. When rotation is turned on, 
the reflection symmetry |(2J) is broken and the pitchforks 
become imperfect, but the SNIPER bifurcation survives, 
with a limit cycle created by one saddle-node bifurcation 
instead of two simultaneous ones. 

Siggers carried out an extensive survey over 4 < 
r < 10, Pr = 0.1 with (less realistic but more tractable) 
stress- free boundary conditions. Like (34|, the numeri- 
cal method used an expansion in eigenfunctions of the 
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FIG. 13: (Color online) Bifurcation diagrams along lines in central (Re, Ra) diagram. T(l/2, 1/2) is plotted as a function of 
Ra for fixed Re in the diagrams along the bottom, and as a function of Re for fixed Ra in the diagrams in the column on the 
right. Diagrams contain steady branches with zero (blue, solid), one (green, long-dashed) and two (red, short-dashed) unstable 
directions. Black dots indicate limit cycles. 
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linearized problem. In [9|, 20-30 eigenfunctions were re- 
tained, a number sufficiently small to allow automatic 
tracking of bifurcations. A phase diagram very similar 
in its large-scale features to that described in this sec- 
tion was found: limit cycles compete with steady states 
that are created by two pitchfork bifurcations and de- 
stroyed by saddle-node bifurcations. However, unlike in 
our case, the basic state never undergoes a Hopf bifur- 
cation like our H; the only Hopf bifurcations are those, 
like our SH, which create small limit cycles, like those in 
inset (i) of figure [12j sometimes called vaccillation. The 
large limit cycle is created by neither a SNIPER nor a 
Hopf bifurcation, but rather by a saddle-node of periodic 
orbits (as in inset (j) of our figure H"2"j) or by a heteroclinic 
bifurcation. 

SNIPER bifurcations are also observed in situations 
with more complicated time dependence, in which the 
saddle-node bifurcations annihilate pairs of limit cycles 
to initiate flow along a torus. This is the case in a 
study of Rayleigh-Benard convection with modulated ro- 
tation [3l| . The basic state is a large-scale axisymmetric 
and reflection-symmetric time-periodic flow. As Ra or 
the modulation amplitude of the rotation is increased, 
pitchfork bifurcations produce pulsating target patterns 
which give way to radially travelling target patterns via 
a SNIPER bifurcation. Abshagen et al. 0, |45| have 
carried out experimental and computational studies of 
secondary bifurcations in Taylor-Couette flow in a small- 
aspect-ratio cylinder. The transitions they consider are 
between nonaxisymmetric rotating and modulated rotat- 
ing waves. However, if one of the time scales is filtered 
out by taking a Poincare map and the reflection sym- 
metry is redefined to include rotation by ir around the 
cylinder axis, their results can be seen as entirely analo- 
gous to ours. Varying the aspect ratio and the Reynolds 
number, the region in which modulated rotating waves 
exist is bounded by curves of SNIPER and Hopf bifurca- 
tions, except over a small region in which homoclinic bi- 
furcations and small-amplitude modulations mediate the 
transition. 



angle 8 while the meridional velocity perturbation is 
v r e r + v z e z = eg/r x V"0. Retaining only the linear 
terms leads to: 

aD 2 ip = -(v r d-+v z d z ) D 2 ^ 

- (K-cL + V z d z ) D 2 i> + (v 2 - D 2 i> 

crG = -(v r d r + v z d z )T 

-(Vrd r + V z d z )e + -^V 2 Q (14b) 
Fr 

av e = - (v r d+ + v z d z ) Vg - (V r d+ + V z d z ) v g 

+ (V - Vg (14c) 

with homogeneous boundary conditions: 

ip = d z ip = <^=> iv = v z = 0, 

vg = 0, 6 = at z = 1 (15a) 
ip = d z ip = <*=^> v r = v z = 0, 

Vg = 0, 9 = at z = (15b) 
ip = d r ip = <^=> v r = v z = 0, 

v g = 0, d r e = at r = 1 (15c) 
ip = D 2 ip = <^=> v r = d r v z = 0, 

vg = 0, d r G = at r = (15d) 

The Reynolds number appears via the inhomogeneous 
boundary conditions on the nonlinear equations (|Si,b), 
making Vg proportional to Re. In order to make this 
dependence on Re explicit, we will scale Vg (but not V r<z 
or *) by Re: 

Vg = Re Ug (16) 



VIII. EIGENVALUES AND EIGENVECTORS 

We now focus on the behavior of the eigenvalues 
and eigenvectors of the basic flow as the Rayleigh and 
Reynolds numbers are varied. Although the linear sta- 
bility analysis of a problem is usually presented before its 
bifurcation scenario, we will focus here on some features 
of the eigenvalues and eigenvectors that are independent 
of the bifurcation scenario. 

To obtain the linear stability equations, we substitute 



(*,V B ,T){r, z) + e^(^,vg,Q)(r, z) 



(13) 



into the governing equations (JS])-©, where (^,Vg,T) 
is a steady solution of ([S])-®, and (ip,vg,@) is a per- 
turbation with growth rate a. We note that the tem- 
perature perturbation is unrelated to the azimuthal 



This scaling can be used even when Re = since Vg 
is proportional to Re. We then rewrite (|14j) in a more 
compact matrix form as follows: 




(17) 



The operators C^, . . . which comprise £, and which de- 
pend on Re implicitly (and weakly) via the steady state, 
are listed below. The boxes surround terms which are 
non-zero when Re = and the steady state is the con- 
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FIG. 14: (Color online) Real part a of leading eigenvalues as a function of Ra for six values of Re. For Re = 0, the eigenvalues 
are real and cross transversely. Zero crossings at Ra — 2260 and Ra = 6640 correspond to pitchfork bifurcations. Purely 
thermal and azimuthal eigenvalues at a — —9.87 and a — —24.6, respectively, are independent of Ra. 

When Re > 0, the transverse crossings become complex conjugate pairs, over Ra intervals which widen with increasing Re. By 
Re = 96, the bifurcating eigenvalue is a complex conjugate pair, leading to a Hopf bifurcation at Ra = 11 856. 



ductive solution V = 0,T = 1 — z. 



-((d z i,)d--(d r ij)d z ) DH 
r 



- (V r d- + V z d z ) £> 2 V> + 



Ra 



-d z (Ugvg) 
r 



-{d z ^)d r T- 
r 



~{d r il))d z T 
r 



--{{d z ^)d + ~{d r ^)d z ) u e 

r 



-( Vr d r + V z d z )G + 
Cv e v e ve = - (V r d + + V z d z ) v g + 



Ft 



D 2 ip 



(18a) 
(18b) 

(18c) 

(18d) 

(18e) 

(18f) 



Figure [TJ] shows the real part of the leading eigenval- 
ues (in units of the inverse viscous diffusive time) as a 
function of Ra for several fixed values of Re, while figure 
1151 shows some of the corresponding leading eigenvectors 
for Re = and for Re = 96, both for Ra = 10 000. Fig- 
ure [T3] shows several eigenvalues which are independent 
of Ra when Re = 0. The origin of these is well known 
and easily understood. For Re = 0, equations (|17p - (fl~8|) 
show that the linear stability problem becomes: 



_I 9r i v 2 
r 1 Er 



/ D 2 







f i> 





I 


i) 


1 8 


V 







\ ve 









( i> 







-) 


e 



(v 



\ Vg 



(18g) with the homogeneous boundary conditions (|15|) . Purely 
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FIG. 15: (Color online) Leading eigenvectors at Ra — 10 000. Left column: Re = 
to convection with one and two toroidal rolls, respectively, with a — 38.5 and a 



and o = -24.6 : 



0. a,b) eigenvectors responsible for transition 
= 15.9. c,d) thermal and azimuthal velocity 



eigenvectors, respectively, with a — —9.87 « 

with a = —21.2. Right column: Re = 96. f,g) first complex conjugate pair, with a ± icu — 
conjugate pair, with a ± i ui = —19.2 ± i 4.35. 



e) eigenvector with two vertically stacked rolls, 
2.24 ± i 3.45. h,i) second complex 



thermal eigenvectors have v = 0, leading to: 
= —Ezg "I 

*e = %v 2 e l=> e = sin (^) ^ 2 

o = eu = a r e| r . 0)1 J = 

(20) 

The first of these thermal eigenvalues is seen in figure [14] 
(i?e = 0) as the horizontal line very close to its analytic 
value of — 7T 2 = —9.8696. The corresponding thermal 
eigenvector is shown in figure [L3b . Another set of solu- 
tions contains only azimuthal velocity: with = tp = 
and vg a solution to 

av e = (V 2 - ^) v e 1 « e = sin(fc7rz) Ji(rj ln ) 

= «eU=o,i = u«|r=o,i J with cr = -((/s7r) 2 + j 2 „) 

(21) 

where Ji is the first Bessel function and j\ n — 3.8317, 
7.0156, ... is one of the zeros of J\. The first of these 



azimuthal velocity eigenvalues is the horizontal line seen 
in figure Q3] (Re = 0), again very close to its analyti- 
cally computed value of -(tt 2 + 3.8317 2 ) = -24.551. The 
corresponding azimuthal velocity eigenvector is shown in 
figure rT5t i- 

The other eigenvalues shown in figure Q3] (Re = 0) 
increase with Ra. The zero crossings of the two largest 
eigenvalues are associated with the pitchfork bifurcations 
discussed extensively in the previous section; these take 
place at Ra — 2260 and Ra = 6640. Their associated 
eigenvectors are shown in figure [L3k.b. for Ra = 10 000, 
where their eigenvalues are 38.5 and 15.9, and they are 
seen to contain one and two concentric radial rolls, re- 
spectively. The leading eigenvectors for higher Re are 
similar. The convective states in figures [5] and [S] do 
not resemble figure 115a . even though they result from 
a pitchfork bifurcation involving this eigenvector. This 
is because the nonlinear steady states in figures [5] and [S] 
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are a superposition of these eigenvectors and the basic 
flows shown in figures |3] and 21 The eigenvector in figure 
[TSb . with two vertically stacked rolls, is associated with 
eigenvalue a = —21.2, between that of the thermal and 
azimuthal eigenvalues at Re = 10 000. 

An overall feature of the eigenvalues for Re = that 
can be seen in figure Q3] is that they are all real and cross 
transversely. Indeed, it is well-known that for Rayleigh- 
Benard convection (i.e. Re = 0), all of the eigenvalues are 
real. When Re > 0, the eigenvalue crossings which occur 
for Re = become complex conjugate pairs. (Because 
the real parts of the eigenvalues are shown in figure [TJ1 
complex conjugate pairs appear as a fusion of two eigen- 
value curves.) As Re increases, the Ra intervals over 
which the eigenvalues are complex widen. By Re = 96, 
the real eigenvalues responsible for the pitchfork bifurca- 
tions have merged into a complex conjugate pair whose 
real part crosses zero at Ra — 11 856 at the Hopf bifur- 
cation point. At Re = 95 (not shown), the fusion of two 
real eigenvalues into a complex conjugate pair and the 
zero-crossing occur simultaneously at Ra = 11750: this 
is the Takens-Bogdanov codimension-two point. 

Figure fTST . g shows the real and imaginary part of the 
complex conjugate pair which is the leading eigenvector 
for Re = 96. The fact that this complex pair at Re = 96 
originates in the fusion of the leading real eigenvectors 
at lower Re is made strikingly clear when the stream- 
functions ip in figure I15f . g are compared with those of 
figure [T5k. b. Note that the choice of the real and imag- 
inary parts of an eigenvector is arbitrary, since an eigen- 
vector can be multiplied by any complex number. The 
particular choice here is imposed by the normalization 
Xm(tp{r = 1/2, z — 1/2)) = 0. The decomposition of 
the complex eigenvector -0 into two vectors containing 
one and two concentric radial rolls shows that the limit 
cycle will involve competing radial structures. As with 
the steady states, the lack of resemblance between figures 
ITST . g and the limit cycle in figure [7] is due to the fact 
that these are superpositions of the basic state and the 
eigenvectors. 

We investigate the progression of the eigenvalues from 
real to complex as Re increases by examining the matrix 
in (fTT)) , which is block diagonal for Re = 0, as shown in 
. Its eigenvalues and eigenvectors thus consist of two 
sets: those of the thermal convection problem (upper left 
submatrix) and those of the azimuthal problem (lower 
right submatrix) . The convective, thermal and azimuthal 
eigenvalues cross transversely because (fT9"]l has no off- 
diagonal terms and because the eigenvalues within each 
of these sub-problems, each associated with a different 
spatial structure, do not cross one another. 

This behavior resembles that seen near the Takens- 
Bogdanov point which occurs in binary fluid or thermoso- 
lutal convection. In the binary/thermosolutal case, the 
pair of eigenvectors which interact have different origins: 
one can be viewed as arising primarily from thermal con- 
vection and the other from solutal convection (46j. In 
the Rayleigh-Benard/von Karman case studied here, the 



eigenvectors which become complex both arise from ther- 
mal convection; the difference between them is their spa- 
tial structure, shown in figure IT5k.b. 

In the binary/thermosolutal case, the transverse cross- 
ings undergo two different fates, depending on the sign of 
the separation parameter 5, which describes whether the 
thermal and solutal convection act in concert or in op- 
position. For positive S, the eigenvalue curves separate 
into two hyperbolas, in what is called avoided crossings, 
and remain real. For negative 5 1 , the eigenvalues join in 
complex conjugate pairs, as in figure [T^l These two cases 
can be understood in terms of a 2 x 2 matrix, whose off- 
diagonal terms are of the same sign if S is positive and 
of opposite signs if S is negative. Here, we have not at- 
tempted such an analysis, but we can conclude that the 
off-diagonal matrices 

( Re ^^ and (ite£,,«, o) (22) 

are in some sense of opposite signs, since the coupling 
they cause between the convective and azimuthal velocity 
eigenvectors leads to complex eigenvalues. The coupling 
between the convective and the purely thermal eigenvec- 
tors must also be of this type. In [46|], a calculation of 
the sign of the coupling terms is presented for the bi- 
nary/thermosolutal case, involving projecting onto the 
eigenvectors of the diagonal submatrices. 

IX. DISCUSSION 

We have shown that a transition analogous to the 
onset of convection occurs in an axisymmetric cylindri- 
cal container subjected to vertical gradients in tempera- 
ture and azimuthal velocity, i.e. in Rayleigh-Benard/von 
Karman flow. Differential rotation, which causes mix- 
ing via Ekman pumping, delays the onset of convec- 
tion. The transition is a pitchfork bifurcation, leading 
to steady convection, for Re < 95 and a Hopf bifur- 
cation, leading to oscillatory convection, for Re > 95. 
Between these two types of convection, over most of the 
(0 < Re < 120,0 < Ra < 30 000) parameter space, the 
transition occurs via a SNIPER bifurcation, in which the 
stable steady states meet a pair of unstable steady states 
and mutually annihilate, leaving a limit cycle in their 
wake. Over a small portion of the parameter space, the 
scenario is more complicated and involves several global 
bifurcations. ( Section II VI mentions some consequences of 
relaxing the imposition of axisymmetry. ) 

The linear stability analysis of the axisymmetric 
Rayleigh-Benard/von Karman problem also shows inter- 
esting features. We have traced the way in which differ- 
ential rotation couples the eigenvalue branches, which are 
real for Rayleigh-Benard convection, in such a way that 
they become complex. The close resemblance between 
the leading real pair of eigenmodes for Re < 95 and the 
leading complex pair for Re > 95 supports the idea that 
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a common basis of eigenvectors could be used to make 
the reduced model studied in the companion paper [34TJ 
more economical. 

Since its discovery by Andronov and Leontovich fill ], 
the SNIPER bifurcation has appeared in a number of 
ODE systems used to model chemical reactions, notably 
in excitable media and the Belousov-Zhabotinsky reac- 
tion [47] . and population biology, for example predator- 
prey systems. In the hydrodynamic and PDE context, 
the SNIPER bifurcation was observed in simulations of 
axisymmetric Rayleigh-Benard convection 0, d, H3, HH . 
Convective states with different radial wavelengths com- 
pete, much as occurs in the present Rayleigh-Benard/von 
Karman configuration. The closest analogy with the 
bifurcation scenario we observe is found in the small- 
aspect-ratio study of Taylor- Couette flow by [441 . l45j . in 
which rotating and modulated rotating waves play the 
role of our steady states and limit cycles, respectively. 

Recently, the SNIPER bifurcation has been the sub- 
ject of renewed attention as a possible explanation for 
the reversals of the earth's magnetic field. A dynamo 
engendered by bulk fluid motion was first produced in a 
laboratory experiment of the von Karman flow in sodium 
(VKS), an electrically conducting fluid This VKS 

experiment shows reversals of the polarity of the mag- 
netic field which bear some similarity to that of the ter- 



restial field. An explanation involving a SNIPER bifur- 
cation in a low-dimensional dynamical system with noise 
has been put forward to explain these reversals [2(| [HJ ■ 
This manifestation of the SNIPER bifurcation also pro- 
vides some possible justifications for studying the ax- 
isymmetric flow. In the VKS experiment, the mean of 
the highly turbulent flow is axisymmetric, even though 
the instantaneous flow is not. A different mechanism 
for magnetic field reversals has been proposed [Hj], in 
which noise is added to a low-dimensional dynamical sys- 
tem displaying a Takens-Bogdanov point. The Takens- 
Bogdanov point, separating a steady from a Hopf bifur- 
cation point, also plays an essential role in our system, 
as it does in the investigations of axisymmetric Rayleigh- 
Benard convection by [9| and of Taylor- Couette flow by 
[45[ • It seems plausible that the TB point and the hys- 
teresis point terminating the SNIPER bifurcation curve 
form part of the unfolding of a codimension-three point. 

Our study was first intended to determine the effect 
of the von Karman flow on the onset of Rayleigh-Benard 
convection. Over the course of the investigation, we en- 
countered a bifurcation scenario leading to oscillations, 
which should and does occur quite generally, whenever 
consecutive pitchfork bifurcations take place. It should 
therefore be of interest to the dynamical-systems com- 
munity as well as to fluid dynamicists. 
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